まず利用するライブラリを読み込みます。Pythonの機械学習ライブラリ sklearn などを利用します。 sklearn については https://scikit-learn.org/stable/index.html を参考にしてください。その他のライブラリについても 勉強してください。
import os # for makedirs
import homcloud.interface as hc # HomCloud
import homcloud.plotly_3d as p3d # 3D visualizaiton
import plotly.graph_objects as go # 3D visualizaiton
import numpy as np # Numerical array library
from tqdm import tqdm_notebook as tqdm # For progressbar
import matplotlib.pyplot as plt # Plotting
import sklearn.linear_model as lm # Machine learning
from sklearn.decomposition import PCA # for PCA
from sklearn.model_selection import train_test_split
%matplotlib inline
%load_ext autoreload
%autoreload 2
まずはデータがどんなものか確認していきます。チュートリアル用のデータは pc というディレクトリ以下にあります。
ls pc/
0000.txt から 0199.txt まで 200個のテキストデータがあります。また、label.txtには 0 か 1 のラベルデータがあり、
これは対象のデータが2つのグループに分けられていることを示しています。
では、ラベルを読み込んで表示します。
labels = np.loadtxt("pc/label.txt")
labels
これは 200個の 0/1の配列です。では、どのデータが0/1でラベル付けされているのかを確認しましょう。 numpy の機能を使います。
np.nonzero(labels == 0), np.nonzero(labels == 1)
0 から 99 までが 0 で、100から199までが 1 で、それぞれラベル付けられていることがわかります。
では、0番目のデータを100番目のデータを可視化してみましょう。homcloud.paraview_interfaceという3次元可視化のモジュールが
HomCloudにはあるので使います。
go.Figure(data=[p3d.PointCloud(np.loadtxt("pc/0000.txt"))], layout=dict(scene=p3d.SimpleScene()))
go.Figure(data=[p3d.PointCloud(np.loadtxt("pc/0100.txt"))], layout=dict(scene=p3d.SimpleScene()))
ではファイルからポイントクラウドデータを読み込みます。
pointclouds = [np.loadtxt("pc/{:04d}.txt".format(i)) for i in tqdm(range(200))]
0番目のデータの「形」を見ると 1000x3 となっていることがわかります。これは3次元上の点が1000個あることを意味します。 どのポイントクラウドも同じです。パーシステント図による機械学習では点の個数は同じくらいでないとあまりうまくいかない (不可能ではないが注意点が多くなる)ことに注意しましょう。
pointclouds[0].shape
では、HomCloudでパーシステント図を計算します。
計算したファイルは pd というディレクトリに保存しておきます。
保存しておいたデータは再利用できるので、別の解析をする場合などにはいちいち計算しなおすことなく利用できます。
os.makedirs("pd", exist_ok=True) # 保存先のディレクトリを作成
for i in tqdm(range(200)):
hc.PDList.from_alpha_filtration(pointclouds[i], save_boundary_map=True,
save_to="pd/{:04d}.pdgm".format(i))
このチュートリアルでは2次のパーシステント図を使います。つまり空隙構造に注目していきます。 もちろん1次や0次なども使えますが、実はこのデータは2次の所に一番特徴があるのでこうします。
より高度な解析法として、1次と2次を同時に使う方法もあります。
pds = [hc.PDList("pd/{:04d}.pdgm".format(i)).dth_diagram(2) for i in tqdm(range(200))]
とりあえず計算したパーシステント図を表示してみます。0番目(ラベル0)と100番目(ラベル1)を表示します。
pds[0].histogram(x_range=(0, 0.03)).plot(colorbar={"type": "log"})
pds[100].histogram(x_range=(0, 0.03)).plot(colorbar={"type": "log"})
Persistence Image については https://arxiv.org/abs/1507.06217 http://www.jmlr.org/papers/volume18/16-337/16-337.pdf を参考にしてください。 ただし、ここで説明するやりかたはこの論文の手法そのままではなく、 https://arxiv.org/abs/1706.10082 https://link.springer.com/article/10.1007%2Fs41468-018-0013-5 で用いられている改変版です。
まずはベクトル化の仕様を決めます。
(0, 0.03), 128 で [0, 0.03]x[0, 0.03] という正方形内を128x128に分割して計算することを指定します。
sigma=0.002 はガウス分布の標準偏差を指定します。また、weight=("atan", 0.01, 3) で対角線からの距離による重みの変化を
指定します。この場合の重みを決めるルールは $\textrm{atan} (0.01 \cdot \ell^3)$ です。ただし$\ell$はlifetime(=birth-death)を意味
します。
vectorize_spec = hc.PIVectorizeSpec((0, 0.03), 128, sigma=0.002, weight=("atan", 0.01, 3))
ではパーシステント図をベクトル化しましょう
pdvects = np.vstack([vectorize_spec.vectorize(pd) for pd in pds])
ベクトルの要素の最大最小を見ると$10^{-10}$と非常に小さいです。これは機械学習に悪い影響を与えるので正規化しておきます。
pdvects.min(), pdvects.max()
pdvects = pdvects / pdvects.max()
まずはPCAを使って解析します。PCAは高次元のデータをできるだけ情報を損わずに低次元に落とすための手法です。 あまり解析結果が破綻しないのでとりあえずパーシステント図との組み合わせでは試す価値があります。 パーシステント図との組み合わせにおいては、個々のパーシステント図を低い次元の点として表現でき、パーシステント図互いの 関係を理解しやすく表現します。
ここでは可視化しやすい2次元に落とします。
pca = PCA(n_components=2)
pca.fit(pdvects)
さてどのように2次元上に表現されているか可視化します。 赤がラベル0、青がラベル1です。
reduced = pca.transform(pdvects) # すべてのデータを2次元に落とす
plt.gca().set_aspect('equal') # 縦横のアスペクト比を揃える
plt.scatter(reduced[labels == 0, 0], reduced[labels == 0, 1], c="r") # ラベル0のデータを赤で描画
plt.scatter(reduced[labels == 1, 0], reduced[labels == 1, 1], c="b") # ラベル1のデータを青で描画
2次元平面上で赤と青がそれなりに分離できているようなので、ここでPIで計算したベクトルは2つのグループの間の違いをうまく表現できているようです。
ではこのX軸とY軸は何を表現しているのでしょうか?ここでは普通のPCAを使っているので、元の次元の所に ベクトルがあり、各点がそのベクトルの成分をどれだけ含んでいるかを表しています。数式で書くと、 X軸、Y軸に対応するベクトル$v_1, v_2$があり、PIによる高次元ベクトル$w_i$($i$はデータごとのインデックスです)は $$ w_i \approx \mu_{i1} v_1 + \mu_{i2} v_2 + v_0 $$ ただし$v_0$は$w_1,\ldots,$の平均、と近似されます(このあたりはPCAの理論を勉強してください)。 すると、$v_1, v_2$がどのようなベクトルであるかがわかればこの図の意味がわかるわけです。
sklearn の PCA では、$v_0$ が PCA.mean_, $v_1, v_2$ が PCA.components_ に保存されています。
pca.mean_
pca.components_
PIの一つの利点として、ベクトルをパーシステント図に変換できる点があります。これを使うとこれら平均や$v_1, v_2$をパーシステント図の形で
可視化できます。実際にやってみましょう。ベクトル化で使ったhc.PIVectorizerMeshのオブジェクトがこの変換(逆変換)のための情報を持っています。
vectorize_spec.histogram_from_vector(pca.mean_).plot()
vectorize_spec.histogram_from_vector(pca.components_[0, :]).plot(colorbar={"type":"linear-midpoint", "midpoint": 0})
vectorize_spec.histogram_from_vector(pca.components_[1, :]).plot(colorbar={"type":"linear-midpoint", "midpoint": 0})
PCAのベクトルの成分には±があり、それが各軸の±と対応しているので±が赤と黒に対応するよう表示させるとわかりやすいです。
それは、colorbar={"type":"linear-midpoint", "midpoint": 0}の部分で調整できます。ここからさらに
この赤い部分、青い部分、に含まれるbirth-death pairを取り出して、逆解析で取り出してさらなる解析をする、ということもできます。
この逆解析の部分は次のロジスティック回帰の部分で実習しましょう。
次はロジスティック回帰です。ロジスティック回帰もPHとの組み合わせが比較的やりやすい手法です。 特に2クラス分類が解釈のしやすさもあってやりやすいです。
まず学習をする前にテスト用データ(test set)と学習用データ(training set)に分けます。
pdvects_train, pdvects_test, labels_train, labels_test = train_test_split(pdvects, labels, test_size=0.25)
学習用データで学習します。sklearnのLogisitcRegressionというクラスを使います。正則化パラメータCはここでは
あらかじめ適当なものを選んでいますが、LogisitcRegressionCVというクラスを使って交差検証を使って自動的にパラメータを
決めることもできます。
model = lm.LogisticRegression(C=0.01, solver="lbfgs")
model.fit(pdvects_train, labels_train, )
テスト用データを使って性能を確認します。
model.score(pdvects_test, labels_test)
なかなか良さそうな精度です。では、学習結果を可視化します。
vectorize_spec.histogram_from_vector(model.coef_.ravel()).plot(colorbar={"type": "linear-midpoint", "midpoint": 0.0})
この図の見方ですが、濃い赤の領域に沢山birth-death pairがあるデータはラベル1である可能性が高く、 濃い青の領域に沢山birth-death pairがあるデータはラベル0である可能性が高いことを意味します。
ではこの赤/青の領域を抽出しましょう。 値が0.005を基準として、0.005より大、-0.005より小、の部分を切り出します。 これには numpy の不等号演算が役立ちます。
coef = model.coef_.ravel()
red_area = vectorize_spec.mask_from_vector(coef > 0.005)
blue_area = vectorize_spec.mask_from_vector(coef < -0.005)
取り出した領域を可視化します。
blue_area.plot(title="blue", ax=plt.subplot(1, 2, 1))
red_area.plot(title="red", ax=plt.subplot(1, 2, 2))
取り出した領域に含まれているbirth-death pairを0番目のPD(ラベル0)100番目のPD(ラベル1)から取り出します。
pairs_in_blue_0 = blue_area.filter_pairs(pds[0].pairs())
pairs_in_red_0 = red_area.filter_pairs(pds[0].pairs())
pairs_in_blue_100 = blue_area.filter_pairs(pds[100].pairs())
pairs_in_red_100 = red_area.filter_pairs(pds[100].pairs())
個数を表示しましょう。
len(pairs_in_blue_0), len(pairs_in_red_0), len(pairs_in_blue_100), len(pairs_in_red_100)
確かに0番目PDには青い領域にあるpairが多く、100番目PDには赤い領域にあるpairが多いようです。実際に見てみると以下のような pairのようです。
pairs_in_blue_0
さて、ここでoptimal volumeを計算して可視化してみます。 詳しくは説明しませんが、
などを参考にしてください。これによってどのような空隙が分類に重要な役割を果たしているのか可視化できます。
optimal_volumes_blue_0 = [pair.optimal_volume(cutoff_radius=0.4) for pair in pairs_in_blue_0]
optimal_volumes_red_100 = [pair.optimal_volume(cutoff_radius=0.4) for pair in pairs_in_red_100]
go.Figure(data=[
ov.to_plotly3d_mesh(color="blue") for ov in optimal_volumes_blue_0
] + [
p3d.PointCloud(pointclouds[0])
], layout=dict(scene=p3d.SimpleScene()))
go.Figure(data=[
ov.to_plotly3d_mesh(color="red") for ov in optimal_volumes_red_100
] + [
p3d.PointCloud(pointclouds[100])
], layout=dict(scene=p3d.SimpleScene()))
ちょっとわかりにくい気もするので2つ重ねて表示してみます。
fig = go.Figure(data=[
ov.to_plotly3d_mesh(color="blue") for ov in optimal_volumes_blue_0
] + [
ov.to_plotly3d_mesh(color="red") for ov in optimal_volumes_red_100
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
青い空隙がラベル0に「典型的」な空隙で、赤い空隙がラベル1に「典型的」な空隙です。
青のほうが赤よりちょっと細い感じがすると思います。学習結果のパーシステント図を見ると、 青い領域は赤い領域よりY軸(death)で小さいほうにあります。実は2次のパーシステント図において、 deathはその空隙に収まる最大の球の半径の2乗に対応しています。つまりこの学習結果は ラベル1のほうがこのスケールの隙間が大きい、ことを意味しています。 赤い領域の中心のdeathが0.022、青い領域の中心が0.017くらいですので、 それぞれの平方根、つまり
np.sqrt(0.017), np.sqrt(0.022)
あたりがそれぞれのグループの特徴的空隙のスケールである、ということがわかります。 元々のパーシステント図を眺めると、もっと小さい空隙はどちらも共通に沢山ある、ということも 見てとれます。
以上でこのチュートリアルは終わりです。